Partly occupied Wannier functions: Construction and applications 
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We have developed a practical scheme to construct partly occupied, maximally localized Wannier 
functions (WFs) for a wide range of systems. We explain and demonstrate how the inclusion 
of selected unoccupied states in the definition of the WFs can improve both their localization 
and symmetry properties. A systematic selection of the relevant unoccupied states is achieved by 
minimizing the spread of the resulting WFs. The method is applied to a silicon cluster, a copper 
crystal and a Cu(100) surface with nitrogen adsorbed. In all cases we demonstrate the existence of 
a set of WFs with particularly good localization and symmetry properties, and we show that this 
set of WFs is characterized by a maximal average localization. 

PACS numbers: 71.15.Ap, 31.15.Ew, 31.15.Ri 



I. INTRODUCTION 

A characteristic property of the single-particle eigen- 
states of most molecular and solid state systems is their 
delocalized nature. For many practical purposes this 
property is undesired and the construction of equivalent 
representations in terms of localized orbitals becomes an 
important issue. 

Within the independent-particle approximation the 
use of Wannier functions (WFs) allows for an exact de- 
scription of the electronic groundstate in terms of a min- 
imal set of localized orbitalsi. The Wannier basis is truly 
minimal in the sense that the number of orbitals is just 
enough to accomodate the valence electrons of the sys- 
tem. Moreover, these localized WFs provide a formal 
justification of the widely used tight-binding^ and Hub- 
bard models 3 . Being the local analogue of the extended 
Bloch states of solid state physics, the WFs formalize 
standard chemical concepts such as bonding, coordina- 
tion and electron lone pairs. Among the more technical 
applications of Wannier functions we mention the con- 
nection to polarization theorji^ and their use within 
so-called "linear scaling" or "order-A^" methods to ob- 
tain the electronic groundstate^. Very recently numerical 
methods for electron transport calculations employing a 
Wannier function basis set have been developed^. 

In the context of molecular systems the analogue of 
Wannier functions for finite systems has been studied un- 
der the name " localized molecular orbitals"^i2iiiii2ii&ii. 
These are traditionally defined by an appropriate unitary 
transformation of the occupied single-particle eigenstates 
and have been used for investigation of chemical bond- 
ing. In the following we shall for simplicity use the term 
WF to cover also localized molecular orbitals. 

In 1997 Marzari and Vanderbilt developed a scheme 
to perform practical calculations of maximally localized 
Wannier functions for an isolated group of bands, i.e. a 
set of bands which is separated by a finite gap from all 
higher- and lower-lying band*^. Within this scheme, the 
usual arbitrariness inherit in the definition of the Wan- 



nier functions due to the unspecified set of unitary trans- 
formations of the Bloch states at every wave vector, is 
removed by requiring that the sum of second moments 
of the resulting WFs is minimal. The method follows 
the traditional idea of defining Wannier functions by a 
unitary transformation of the occupied (Bloch-) orbitals. 
In general, such methods fail to produce well localized 
orbitals when applied to metallic systems because the 
unoccupied states belonging to the partly filled valence 
bandsi& are not considered. Of course, in cases where the 
partly filled valence bands are separated by a gap from 
all higher bands, the method of Marzari and Vanderbilt 
still applies. However, in the more general case where 
the bands of interest cross and/or hybridize with other 
unwanted bands a different approach must be used. 

In this paper we demonstrate how the localization and 
in some cases also the symmetry of a set of WFs can 
be drastically improved by including selected unoccupied 
states in the definition of the WFsi&. The determina- 
tion of the relevant unoccupied states can be viewed as a 
bonding-antibonding closing procedure, where occupied 
bonding states are paired with their antibonding coun- 
terparts to yield localized orbitals. To be more specific, 
consider two well-localized atomic orbitals on neighbor- 
ing atoms in a molecule. If we allow the two states to hy- 
bridize, a bonding and an antibonding combination will 
result - combinations which may be less localized than 
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FIG. 1: Schematic of the bonding-antibonding closure for a 
hydrogen molecule. The construction of well-localized atomic 
s-orbitals involves a matching of bonding and antibonding 
orbitals, independent of their occupation. The sign of the 
wave functions is indicated by the shading. 
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the individual atomic orbitals. To regain the localized 
atomic orbitals from the molecular orbitals we need both 
the bonding and the antibonding combination indepen- 
dent of their occupation, see Fig. ^ In some cases the 
antibonding state may have hybridized further with other 
states and the state which "matches" the bonding state 
will be a linear combination of eigenstates. The prob- 
lem we address here is the construction of a method for 
systematically identifying the relevant unoccupied states. 
We show that this can be achieved by optimizing the lo- 
calization of the resulting WFs. The paper gives a more 
detailed and extended account of the work previously 
published in a Letter 

For periodic systems the bonding-antibonding closure 
can be viewed as a procedure for disentangling the partly 
occupied valence bands from higher-lying bands. This 
problem has previously been addressed by Souza et alw& 
who proposed a disentangling method based on a mini- 
mization of the change in character of the Bloch states 
across the Brillouin zone (BZ). While this is a natural 
strategy for crystalline systems, it is not clear how this 
disentanglement procedure applies to non-periodic sys- 
tems like isolated molecules, a surface with adsorbates or 
a metal with impurities. 

The present method is related to that of Souza et ali&, 
however, instead of minimizing the dispersion across the 
BZ we suggest a disentanglement procedure based exclu- 
sively on a minimization of the spread of the WFs. In 
this way we omit any reference to the wave- vector and are 
therefore not limited to periodic systems. The generality 
of the method is demonstrated by application to three 
different systems: an isolated Sis cluster, a copper crys- 
tal, and a Cu(100) surface with nitrogen adsorbed. Our 
results for the copper crystal are very similar to those ob- 
tained by Souza et alm&, and this indicates the similarity 
of the two localization schemes for periodic systems. 

The paper is organized as follows: In Sec. ^ we in- 
troduce the spread functional and outline the strategy 
behind the localization algorithm. In Scc. lIIII we give the 
formal definition of partly occupied WFs in the limiting 
case of a large supercell and derive the corresponding ex- 
pressions for the gradient of the spread functional. The 
extension to periodic systems is discussed in Sec. IIVI In 
Sec. [V] we apply the method to a Sis cluster, a copper 
crystal and a Cu(100) surface with adsorbed nitrogen. 



II. DESCRIPTION OF THE METHOD 

In this section we introduce the spread functional used 
to measure the degree of localization of a set of orbitals, 
and give an introductory description of the localization 
scheme including its relation to the method of Souza et 
alii 



A. Spread functional 

Within the localization scheme of Marzari and Van- 
derbilt 1 ^ the spread of a set of functions {w n (r)}^ =1 is 
measured by the sum of second moments 

N 

s = ]T(K|r 2 K> - K|r| W „) 2 ). (1) 

71=1 

When periodic boundary conditions are applied, as in the 
present study, and the supercell is sufficiently large, the 
minimization of S is equivalent to the maximization of 21 

N N G 

n = ^£vF Q |Z Q , n „| 2 , ( 2 ) 

n— 1 a — 1 

where the matrix Z a is defined as 

Z a , nm = {w n \e~ lGa ' r \w m ). (3) 

The {G Q } is a set of at most six reciprocal lattice vectors 
and {W a } are corresponding weights which account for 
the shape of the unit cell. For a defi nition and discussion 
of these quantities we refer to Rcfs. 11.31141 

B. The localization scheme 

The starting point is the set of single-particle eigen- 
states, {ipn\, resulting from a conventional electronic 
structure calculation. For simplicity we shall assume that 
the system is isolated or is contained in a large supercell 
such that reference to k-points can be omitted. The aim 
is to obtain a set of N w localized WFs with the prop- 
erty that any eigenstate below a specified energy, E$, 
can be exactly reproduced as a linear combination of the 
WFs. An obvious way to achieve this would be to ap- 
ply the method of Marzari and Vanderbilt to compute 
the unitary transformation of the N w lowest eigenstates 
leading to the most localized WFs. The problem with 
this strategy is, however, that it is in general not possi- 
ble to localize all WFs simultaneously, and the problem 
cannot be overcome by increasing N w . 

Instead, we define an external localization space as 
the space spanned by the Nb lowest-lying eigenstates 
(Ni, > N w ). Within this space we consider the subspace 
spanned by the eigenstates with energy below E , to- 
gether with L extra degrees of freedom (EDF). We shall 
refer to this subspace as the active localization space or 
simply the localization space. The EDF are assumed to 
be orthogonal and L is chosen such that the dimension 
of the active localization space equals N w . We then per- 
form a simultaneous optimization of the WFs within the 
active localization space and of the active localization 
space itself. In practice this is achieved by optimizing an 
N w x N w unitary matrix together with the coordinates of 
the EDF such that the functional il becomes maximal. 
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It is the determination of the EDF that distinguish 
our method from that of Souza et al. 19 In the latter, 
the spread functional is decomposed into two terms: 
n = fti + il, where il/ is related to the fc-space disper- 
sion of the band-projection operator, see Ref. In the 
first step, the EDF are determined by maximizing ilj, 
which depends only on the localization space itself and 
not on the internal unitary transformation. In the sec- 
ond step f2, or equivalently O, is then maximized within 
the fixed localization space. It is clear, that the sepa- 
rate maximization of f2j and f2 does not amount to the 
global maximization of f2 that we propose here. We shall, 
however, see that the two methods lead to very similar 
results in the case of periodic systems. 



III. LARGE SUPERCELLS 

In this section we give a detailed description of the 
localization scheme in the limiting case of a large super- 
cell where a T-point sampling of the first Brillouin zone 
is a good approximation. For simplicity we discuss this 
case separately before extending it to periodic systems, 
although the latter contains the former as a special case. 
After giving the definition of partly occupied Wannicr 
functions we derive expressions for the gradients of the 
spread functional and discuss how to combine these with 
a Lagrange multiplier scheme to determine the maximum 
of £1 



A. Definition of partly occupied Wannier functions 

We denote the total number of eigenstates obtained 
from the electronic structure calculation by Nf, and the 
number of eigenstates below the energy Eq by M . Our 
aim is to construct a set of N w WFs which span at 
least the M lowest-lying eigenstates. The remaining 
L = N w — M degrees of freedom are simply used to im- 
prove the localization of the resulting WFs as much as 
possible. We expand the WFs in terms of the M lowest 
lying eigenstates and L extra degrees of freedom, {<fii}, 
belonging to the (N — M)-dimensional space of eigen- 
states with energy above Eq: 



M 



2J Umn^m + ^ U M 



+l,n<Pl, 



(4) 



1=1 



where the extra degrees of freedom (EDF) are written as 



N b -M 
k = 2J c ml1pM+T. 



(5) 



The columns of the matrix c are orthonormal and 
represent the coordinates of the EDF with respect 
to the eigenstates lying above Eq. The matrix U 



is unitary and represents a rotation of the functions 

{ipi,...,ip M ,<i>i,---,<t>L}- 

In order to simplify the notation we introduce the ma- 
trices 



C 



jMxM q 

c 



V = cu 



U M 
cU L 



(6) 



where U M and Ul denotes the M uppermost and L low- 
ermost rows of U 7 respectively. The ith column of V 
gives the coordinates of Wi with respect to the full set of 
eigenstates {ifj n }- 

Substituting the expansions Q and JSJ into Eq. © 
we obtain a compact matrix expression 

Z a = V t z£ V = [ft & CU, (7) 

where is obtained from Eq. J3J by using the eigen- 
states {rjj n } in the inner product, 

4°L=(^|e- iG 



IVv, 



(8) 



Gradient of 



Through Eq. Q the spread functional, tt, in Eq. @ 
becomes a function of the matrices U and c. The max- 
imum of n can be found iteratively by updating U and 
c in the direction given by the gradient. In the following 
we derive expressions for the gradient of f2. 

We write the unitary matrix at iteration n as 
[/("— 1) exp(— A), where A is an anti-hermitian matrix. 
Since we are only concerned with small variations, we 
expand the exponential to first order, i.e. exp(— A) i£ 
1 — A. Inserting this into Eqs. (Q and J5J we find 



on 

dA~ 



N G 



' Z a ij(Za,ii Z a jj 



)]• 



(9) 

All matrices in this expression refer to iteration n — 1. 
The new rotation at iteration n is then obtained by mul- 
tiplying [/( n_1 ) by exp[— c^V^O)] where d is the length 
of the steepest-ascent step and [Va^j = dtt/dAij. 

We now turn to the problem of determining the steep- 
est uphill direction of f2 with respect to variations in c. 
In general, for a real-valued function j(z = x + iy) the 
direction of steepest ascent with respect to z is given by 

9/ 

dz* 



2 [ dx +t dy ) - 



(10) 



To calculate the gradient dft / <9c* ? - we use that 



d\Z a , 



dc* 



— Z a 



dZ* 



dc* 



Z* 



dZ a , 



From Eq. JJJ) it follows that 



dc* 



dZ arLi 
dc* 



U na a * ^ a ,bc cd dn 



na g * 

abed y 



(11) 



(12) 



E 

abed 



r/t r+ z(°) dCcd TJj 

u na'~'ab' a,bc Q _* u dn ' 



dc*- 
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and from definition (jfj) 



= 



5C t 

nm <r r 



It is now easy to establish that 



dc* 



[z^ ) v] M + l , n u* M+:jn 



dZ* 



dc* 



[(ZM)W] M+ i,nU 



M+j,n- 



(13) 
(14) 

(15) 
(16) 



Combining Eq. I|ll|) with l|15|) and l|ltj|) we arrive at the 
desired expression 



on 

dc*~ 



Y^W a [zWVB(Z* a )U^+(Z^VB(Z a )U^] 

a=l 



M+i,M+jy 



. (17) 
in the 



where B(Z a ) is a diagonal matrix with [Z a< 
diagonal. 

To treat the constraint that the EDF {(j>i} should be 
orthonormal during the maximization procedure we in- 
troduce the Lagrange multipliers Ay- and perform an un- 
constrained maximization of the functional 



(18) 



The Lagrange multipliers are initially unknown and must 
be estimated at each iteration. At the maximum we have 
V c *f2z, = which is equivalent to the condition 



V c .fi- c A = 0. 
Multiplying by cr from the left leads to 
A T = c + V c »0. 



(19) 
(20) 



This relation can be used to estimate the Lagrange multi- 
pliers at each iteration. A step of length d in the steepest 
uphill direction is thus accomplished by adding to c the 
matrix d(l — cc^)V c *£l, followed by an orthonormaliza- 
tion of the columns of c. 



IV. PERIODIC SYSTEMS 

We consider a periodic system with a unit cell defined 
by basis vectors ai , &2 , a3 which in turn define the ba- 
sis vectors of the reciprocal lattice bi, b2, b3. The Bloch 
states, {VVik}, resulting from the electronic structure cal- 
culation are characterized by a band index n and a crystal 
momentum k. The total number of bands is denoted by 
Nb and the number of eigenstates at a given k-point with 
energy below E is denoted by Mk. We assume a uni- 
form sampling of the first BZ such that any k-point can 
be written as 



Ni 



bi 



n 2 
N 2 



n 3 
N 3 



(21) 



where N is the number of k-points in the direction b^ 
and rii = 0, . . . , JYj — 1. Note that the T point is al- 
ways included. With this convention the Bloch states, 
{V'nk}, correspond exactly to the T-point eigenstates of 
the repeated cell defined by the extended basis vectors 
iViai , A^a2 , A?3a3 . An alternative way of stating this cor- 
respondence is to say that the k-points in Eq. (|21|l fall 
on the reciprocal lattice of the repeated cell, see Fig. [21 
As we shall see below, this correspondence allows us to 
use the spread functional fl defined in Eq. @ also for the 
periodic system. We stress that the formalism developed 
in the following section contains the T-point formalism 
described in the preceding sections as a special case. 



A. Definition of partly occupied Wannier functions 

We write the nth Wannier function related to unit cell 

i as 



w, 



(22) 



where N/. is the total number of k-points and ipnk is a 
generalized Bloch state to be defined belowi^. Each gen- 
eralized band, i.e. each set {ipnk} for fixed n, gives rise 
to one WF per unit cell. These WFs are simply related 
by translation, i.e. Wi^ n (r) — wo in (r — R-i), and thus it 
suffices to consider the WFs of the cell at the origin. In 
doing this we can omit the cell index and simply denote 
the WFs by {w n }. We denote the number of WFs per 
cell by N w . 

Following the idea behind Eq. we expand the gen- 
eralized Bloch state t^nk in terms of the Mk lowest lying 
Bloch states and extra degrees of freedom, {0/k}, from 
the remaining (Nb — A/k)-dimcnsional space 

4>nk = ^L^mk + U M k +l,n^> ( 23 ) 

m=l 1 = 1 

where the EDF are expanded as 



N b -M k 

E 

m— 1 



(24) 



The number of EDF at a given k-point is determined by 
the condition Lk + Mk = N w . If Mk exceeds N w , we 
simply put Mk = N w . Due to the exact correspondence 
between the Bloch states {ifink} and the T-point eigen- 
states of the repeated cell, we can use the functional © 
to measure the spread of the Wannier functions. The 
matrices Z a are still defined by Eq. © but it should be 
remembered that the inner product as well as the recip- 
rocal lattice vector G Q now refer to the repeated cell. 
From Eqs. I|23I24|) we find the following generalization of 
Eq. CD 



Za — Y2 



kk' 



(25) 



k,k' 
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where 

Z kk ' = (f/ k ) t (C k ) t Z( t ^ kk 'C k 'C/ k '. (26) 

The matrix C k is given by the obvious k-point analogue 
of Eq. © and the matrix zi°^ kk is defined by 

= (V>«k|e- iG - r |Vw>. (27) 

Most of the matrices zi°' ) ' kk are in fact zero. Writing 
the Bloch functions as i/'nk — %k(r) exp(ik • r), where 
u n k has the periodicity of the lattice, we get 

4%™ ' = J < k (r)« mk ,(r)e l ( k '- k - G °)- r dr, (28) 

which is non-zero only when 

k' = k + G a . (29) 

Here it is implicit that k and k' belong to the first BZ 
and thus it might be necessary to translate k' by a re- 
ciprocal lattice vector. The relation between k and k' is 
illustrated in Fig. Note that the condition in Eq. |2l^l 
reduces the double sum in Eq. I|25|) to a single sum over 
k. 



on 



on 



We note that these expressions, of course, reduce to 
Eqs. I|9I17|I in the limit of a single k-point. The max- 
imization of £7 proceeds along the same lines as for the 
T-point case, except that Lagrange multipliers are needed 
for each k-point. For example the analogue of Eq. I|18|) 
reads 

= & - y^Afj ]k (0ik|^ik)- (32) 

B. Optimizing the number of extra degrees of 
freedom 

For given values of Nf,, N w and Eo, the algorithm in- 
troduced above produces the N w most localized WFs 
that can be formed within the external localization space 
when all eigenstates below E should be exactly repro- 
ducible in terms of the WFs. It remains to determine the 




FIG. 2: Relation between the first BZ of the unit cell, denned 
by the reciprocal basis vectors bi, D2, D3 (light gray), and the 
first BZ of the repeated unit cell (dark gray). In this case iVi 
and N2 from Eq. 12H both equals 3. The relation between k 
and k', given in Eq. I|29|l . is indicated. 



The derivation of the gradient of £1 follows closely the 
T-point case discussed in Sec. IIII El and is therefore omit- 
ted. The result is 



(30) 



(31) 

I 

optimal values for Nf, and N w for a given Eq . Let us start 
by considering the situation where Nb has been fixed at 
a value which is large enough to include all anti-bonding 
states relevant for the localization. In practice this typi- 
cally means ~ 10 eV above the Fermi level. It seems as a 
natural strategy to choose N w such that the localization 
per orbital is maximal. To quantify this condition we 
define the average localization per orbital as 

(O) = (33) 

where we have indicated the dependence of £1 on the 
three parameters explicitly. We note that since the value 
of £1 also depends on the size and shape of the super- 
cell, it does not make sense to compare the value of £1 
for systems described in different supercells. Fixing N w 
on the basis of (£1) represents a completely general cri- 
terion which can be applied in any situation. However, 



£ 

0=1 



T/T7 U 7 H 7K- G a .k . ry /yk,k+G Q \* I ry X^^k.K+G^ ry //7k— G a ,k\*l 



No 



J2 w^a[4 0) ' k ' k+G ^ k+G °D(z:)([/ k )t + (z(°>< k - G - k )ty k - G «D(z Q )([/ k )t; 



M k +i,M k +j- 
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the localization procedure must be carried out for several 
values of N w which might be a tedious task depending 
on the size of the system. We next consider the situa- 
tion when Nb is also allowed to change. Formally, the 
global maximum of (f2) is attained in the limit where 
both Nb and N w tend to infinity in which case an in- 
finite set of completely localized delta functions can be 
realized. However, we have found that for practical val- 
ues of Nb where very high energy states are not included 
in the external localization space, (ft) will have a local 
maximum for some N w , and the position of the maxi- 
mum is not sensitive to the actual value of Nb. Thus, it 
is indeed possible to determine an optimal value of N w 
by maximizing (17). 

Alternatively it is often possible to determine a value 
for N w based on symmetry arguments, chemical intu- 
ition, or a closed band condition. As we shall see in the 
following examples the two criteria for determining N w 
lead to similar results. 

C. Start guess for U k and c k 

For small systems we have found that the localization 
algorithm is quite stable and usually leads to the global 
maximum independently of the initial value of the ma- 
trices {C/ k }, {c k }. For larger systems, however, there is a 
risk of getting stuck in a local maximum and in such cases 
the start guess becomes important. It is then natural to 
start from a set of simple orbitals located either at the 
atoms or at the bond centers. Let {f v } denote such a set 
of simple orbitals. The question is how to transform this 
into the matrices {C^ k }, {c k }. To this end we project the 
initial orbitals onto the subspace spanned by the Bloch 
states at each k-point: 

fuk = / l{lpnk\fu)lpnk- (34) 

71=1 

The following procedure is carried out for each k-point 
separately. For fixed k we regard (ipnk\fv) as a matrix in 
the indices n, v. Its columns represent the coordinates of 
the /kt/ with respect to the Bloch states {^>nk} n =i an d 
as such it is a (non-orthogonalized) version of the matrix 
y k , see Eq. JBJ. After a normalization of the columns of 

(V'nkl/i/) we compute the norm of the component of /„k 
orthogonal to the occupied subspace: 

ll/Jkl| 2 = E K<M/„)| 2 . ( 35 ) 
n=Af(k) 

The first EDF is chosen as a normalized version of the 
/ki/ f° r which H/k^ll is the largest. The remaining /^'s 
are then orthogonalized onto this vector and the pro- 
cess is repeated until all EDF, and thus c k , have been 
determined. Finally the identity C/ k = (C k )^V k with 
— > (ipnk.\fu) determines L/ k . Since the f v u are not 
necessarily orthogonal, the columns of the resulting J7 k 
must be explicitly orthogonalized. 



V. RESULTS 

In the following sections we apply the localization 
scheme to three different systems. To demonstrate the 
generality of the method we consider both isolated and 
metallic systems as well as a metal surface with adsorbed 
impurities. In Sec. IV Al we construct partly occupied 
WFs for an isolated Sis cluster and illustrate how differ- 
ent sets of WFs can be obtained by varying the number of 
extra degrees of freedom. In Sec. IV Bl we investigate the 
WFs of a Cu(fcc) crystal and compare the results with 
those obtained by Souza and co-workers^ who studied 
the same system using a different but related method. 
Finally, in Sec. IVCl we perform a detailed WF analysis 
for a Cu(100) surface with 0.5 mono-layers of nitrogen. In 
all calculations we use a plane- wave based DFT code2£ to 
obtain the Kohn-Sham eigenstates, and we describe the 
ion potential by Vanderbilt ultrasoft pseudopotentials 24 . 
To ensure a proper convergence of the unoccupied states 
employed in the localization scheme, the DFT calcula- 
tions have been converged with respect to the full set of 
Kohn-Sham eigenvalues. In the appendix we explain how 
to extend the localization scheme to ultrasoft pseudopo- 
tentials. 



A. Sis cluster 

As an example of an isolated system we consider an Sis 
cluster in its ground-state geometry^, see Fig. |3f a) . We 
use a cubic supercell of length 16 A and sample the first 
BZ at the T-point. To test the dependence on the size of 
the external localization space we consider the two cases 
N b = 30 and N b = 100. We set M = 10 corresponding to 
the number of occupied states, and calculate the average 
localization per WF, (fi) = Q/N w , for L = 0, . . . , 7. The 
result is shown in Fig.0] For L = there is no difference 
between the two cases since the WFs are constructed en- 
tirely from the occupied eigenstates. However, for L > 
the larger space available for the extra degrees of freedom 
leads to an improved localization when Nb — 100. Apart 
from this general improvement in localization, there is 
no qualitative difference between the WFs obtained with 
N b = 30 and N b = 100 for a given L. We note that both 
curves have a maximum for L — 4, corresponding to a to- 
tal of 14 WFs. This particular set of WFs together with 
their centers is shown in Fig.[3fb). The fact that this set 
of WFs respect the symmetry of the cluster is a special 
property of the L = 4 solution: for other values of L, 
including L — 0, the WFs break the symmetry of the Sis 
cluster. This indicates that the solution corresponding 
to the maximal value of (Q) has a special meaning. In- 
deed, the value N w — 14 could also have been anticipated 
from physical arguments. Starting from a set of four sp 3 
orbitals located at each Si atom we expect bonding and 
anti-bonding states to form between pairs of aligned or- 
bitals belonging to nearest neighbor pairs of Si atoms. 
These bonding states can be identified as the six bond- 
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centered WFs shown to the far left in Fig. [3] The two 
"top" Si atoms have three nearest neighbors and thus a 
single sp 3 orbital is left as a lone pair (middle WF). The 
remaining three Si atoms each have two nearest neighbors 
and consequently two sp 3 orbitals are left as lone pairs 
(rightmost WF). In total this adds up to 14 orbitals. The 
anti-bonding counterparts of the bonding states formed 
between nearest neighbors are not brought into play for 
L = 4, because they are much less localized than the 
bonding states. However, by setting L — 10 and thus 
searching for a total of 20 WFs, the anti-bonding states 
are picked out as EDF and we obtain a full set of sp 3 
orbitals. This solution has, however, a smaller value for 
(fi) than the solution at L = 4. 




x6 x2 x6 



FIG. 3: (a) Geometry of the Sis cluster, (b) Contour plots 
of the WFs corresponding to L = 4. The position of the WF 
centers are indicated by black spheres. 



B. Copper crystal 

To illustrate the method in the case of a periodic sys- 
tem we consider the construction of WFs for a copper 
crystal. This system was also studied by Souza et alm& 
using their disentangling method to obtain the WFs. Our 
results are in close agreement with those obtained by 
Souza et ai, and this indicates the similarity of the two 
methods for periodic systems. 

We use the primitive fee unit cell and sample the first 
BZ on a uniform (11,11,11) Monckhorst pack grid con- 
taining the r-point. To obtain a minimal set of WFs 
describing the Cu d- and s-bands we set N w — 6. We 
construct two sets of WFs corresponding to two differ- 
ent values of E : (i) E = 0.0 eV and (ii) E Q = 3.0 eV, 




1 2 3 4 5 6 7 
Extra degrees of freedom (L) 



FIG. 4: Average spread of the WFs of the Sis cluster for 
different values of Nb and L. 

relative to the Fermi level. In the first case the resulting 
WFs will span at least the occupied subspace and thus 
the electronic structure described by the WFs will be cor- 
rect below Ep. In the second case the electronic struc- 
ture will be correct up to 3 eV above Ep, however, since 
this is a stronger restriction on the localization space we 
must expect that the resulting WFs will be less local- 
ized than those obtained in (i). In Figs. [S] and we 
show the original DFT bands together with the approxi- 
mate bands computed by diagonalizing the Hamiltonian 
within the subspace spanned by the WFs of case (i) and 
(ii) , respectively. In both cases we see a very good agree- 
ment between the exact and approximate bands below 
Eq. At higher energies the approximate bands deviate 
from the exact bands, indicating that the EDF which op- 
timize the localization of the WFs do not coincide with 
specific Bloch eigenstates. The quality of the WF bands 
below Eq depends on the number of k-points used to 
construct the WFs. This is because the band diagram 
must be constructed starting from fully localized func- 
tions, which means that the coupling matrix elements 
must be truncated beyond a cut-off distance given ap- 
proximately by Ni/2 unit cells in the direction a^. Thus 
the repeated cell, or equivalently the number of k-points, 
must be so large that the WFs have decayed sufficiently 
between the repeated images. 

Both sets of WFs consist of five atom-centered d- 
orbitals and a single s-like orbital centered in one of the 
two interstitial sites. The c?-orbitals are more or less iden- 
tical for the two cases, and two examples arc shown in 
Fig. Contour plots of the s-like orbital is shown in 
Fig. [S] (b) and (c) for case (i) and (ii) , respectively. The 
plots indicate that the s-orbital of case (ii) is less local- 
ized than the one obtained in case (i). That this is indeed 
correct follows from the value of the spread functional, 
O, which is higher for (i) than for (ii). 

The minimal set of WFs obtained with N w = 6 breaks 
the symmetry of the fee crystal because the s-like orbital 
is located in one of the interstitial sites leaving the other 




FIG. 5: Band structure of Cu(fcc). The full lines are the orig- 
inal DFT bands and the dots are the approximate bands com- 
puted from a set of six WFs (N w = 6). The WFs have been 
constructed using llxllxll k-points and keeping all occupied 
states in the localization space, i.e. Eo = 0.0 eV relative to 
the Fermi level. 
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FIG. 6: Like Fig.|3except that the WFs have been generated 
with E = 3.0 eV. 




FIG. 7: Two d-like WFs for Cu(fcc). The orbitals are cen- 
tered at the atoms (not shown) . 



FIG. 8: (a) Two tetrahedral interstitial sites in the fee crystal, 
(b-d) Contour plots of the s-like WF obtained for Cu(fcc). 
The WFs in (b) and (c) have been generated with N w — 6, 
E = 0.0 eV and N w = 6, E = 3.0 eV, respectively. Both 
WFs are located in one of the interstitial sites. The WF in 
(d) correspond to N w = 7 and Eo = 0.0 eV. In this case there 
is an equivalent WF located in the other interstitial site. The 
same contour value has been used for all plots. 
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FIG. 9: Like Fig. [HJexcept that the WFs have been generated 
with N w = 7. 



empty. As demonstrated by Souza et al. the symmetry 
can be restored by using seven WFs per primitive cell 
instead of six. In Fig. we show the band structure ob- 
tained from a set of WFs generated with N w = 7 and 
Eq = 0.0 eV. We note that very high-energetic states 
are now selected as the optimal EDF. This solution can 
therefore only be obtained for rather large external lo- 
calization spaces, i.e. Nj, > 9. The five d-like WFs are 
unchanged, but now we obtain two equivalent s-like WFs 
located in each of the two interstitial sites thereby restor- 
ing the fee symmetry, see Fig.|SJd). We have calculated 
the average localization (fi) for N w = 6, 7, 8, and found 
that the maximum is attained for the symmetric solution 



with N w = 7. 



C. Nitrogen absorption on Cu(100) 

In this section we study the WFs of a copper (100)- 
surface covered with half a mono-layer of nitrogen atoms. 
As the system is neither periodic (in all directions) nor 
isolated, it represents a very general situation. The sec- 
tion is divided into two parts. In the first part the WFs 
are constructed and analyzed, and in the second part we 
use the obtained WFs to study the chemisorption of ni- 
trogen within the Newns Anderson model. 
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1. Wannier function analysis 

We model the Cu(100) surface by a slab with a thick- 
ness of two atomic layers. The supercell contains four Cu 
atoms and a single N atom adsorbed in a hollow site, and 
its height is such that the surface slabs are separated by 
9.0 A of vacuum. A topographic top-view of the surface 
is shown in Fig. llUI We sample the first BZ on a uniform 
(7,7,1) Monckhorst pack grid containing the T-point. 




FIG. 10: Topographic view on the Cu(100) surface with ad- 
sorbed nitrogen. The white spheres are nitrogen, while the 
light gray spheres represent the Cu surface layer. A supercell 
is indicated. 




Number of Wannier functions (N ) 



FIG. 11: Average localization of the WFs of the nitrogen 
covered Cu surface for different values of Nb and N w . 

Let us start by considering what we can expect to find 
on the basis of our previous experience. First, the result 
from the copper crystal suggests that a minimal descrip- 
tion of the metal surface is obtained with five d-orbitals 
and an s-like orbital per Cu atom. Since there are four Cu 
atoms per supercell this gives a total of 24 WFs. Next, 
the similarity between the valency of N and Si together 
with our experience from the Sis cluster points to a de- 
scription of the nitrogen atom in terms of sp^-hy brides. 



In Fig.^Jwe have plotted the average localization, (O) , 
of the obtained WFs as a function of N w for three differ- 
ent sizes of the external localization space corresponding 
to Nf, — 35, 40, 50. In all cases we have set Eq — Ep in 
order to ensure that the occupied eigenstates are exactly 
reproduced by the WFs. As expected, the localization 
improves as the size of the external localization space in- 
creases. In addition, the maximum of (Q) shifts towards 
larger L-values as Nb is increased. Specifically the maxi- 
mum shifts from N w — 27 to N w = 29 as Nj, is increased 
from 35 to 50. This is not unexpected since we know 
that (f2) will be a monotonically increasing function of L 
in the limit Nb — > oo, see discussion in Sec. II V Bl Again 
we stress that it is only the degree of localization of the 
WFs that change with Nb for a fixed L, and not their 
qualitative form. Thus the chemical picture provided by 
the WFs does not change with Nb. In fact, for all the 
values of Nj, we obtain 20 highly localized d-orbitals (five 
located on each of the four Cu atoms) and four sp 3 or- 
bitals centered on the N atom, see Fig.^H The remaining 
N w — 24 WFs are the less localized s-like orbitals of Cu. 
Thus, as N w is increased beyond 24, the number of s-like 
Cu WFs simply increases correspondingly. 




FIG. 12: Average localization of the d-, sp 3 - and s-like WFs 
considered separately for different values of Nb and N w . 

To gain further insight into the dependence of the WFs 
on Nb and N w , we show in Fig.^|the average localization 
of the d-, sp 3 - and s-orbitals, separately. It is clear that 
the -/^-dependence as well as the maximum of (f2) are 
almost exclusively related to the Cu s-orbitals. Except 
for the case Nb = 50, which is in fact somewhat extreme 
since states of 20 eV above the Fermi level are included 
in the external localization space, the average spread of 
the Cu s-orbitals is maximal for N w = 28. This corre- 
sponds to one s-orbital per Cu which is exactly what we 
anticipated from the analysis of the copper crystal. 

We end by summarizing the chemical picture obtained 
from the WF analysis: For N w = 28 the Cu surface is 
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described by the minimal set of WFs consisting of five 
d- and one s-like orbital per atom. For the nitrogen we 
obtain four sp 3 hybrids oriented as indicated in Fig. EH 



p° g denotes the PDOS of the group orbital in the absence 
of coupling to the adsorbate state. Often p° g {e) is referred 
to as the band to which the adsorbate is coupled. 




(Down type) 

FIG. 13: Orientation of the four sp 3 -like WFs belonging to 
the N atom. Seen from above the orbitals point to the bridge 
sites of the Cu atoms in the surface. Each orbital points either 
up from or down into the surface. An example of each type 
is shown in the lower panel. 



2. Adsorption in the Newns Anderson model 

The WFs can be used to obtain a detailed and consis- 
tent picture of the hybridization occurring between the 
nitrogen states and the states of the substrate. As we 
shall see the analysis gives a complete account for the 
shape of the projected density of states of a given N 
orbital, in terms of the bare orbital energy, a coupling 
strength, and the density of states of the so-called group 
orbital. 

In the Newns Anderson model, one considers an ad- 
sorbate state, \a), of energy e a = (a\H\a), coupled 
to a continuum of states, \k), representing the sub- 
strate. The coupling matrix elements are denoted by 
Vk = (a\H\k). A particularly useful formulation can be 
obtained by introducing the normalized group orbital, 
Is) = ^EfcW, where V = (E fe l^| 2 )" 1/2 - It is 
easily checked, that the coupling between \a) and any 
substrate state orthogonal to \g) vanishes. Consequently 
| a) is coupled to the substrate via the group orbital only, 
and the coupling is given by V, i.e. V = (a\H\g). 

Physical quantities such as the projected density of 
states (PDOS) of the adsorbate state and the hybridiza- 
tion part of the adsorption energy, can be obtained from 
the retarded adsorbate Green's function, which in turn 
follows from the three quantities e a , V and p„{e), where 




-4 -2 
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FIG. 14: Top: PDOS for the atomic p-orbitals obtained by 
diagonalizing the Hamiltonian in the subspace spanned by 
the sp 3 WFs of the N atom. Bottom: PDOS for the group 
orbitals corresponding to each of the atomic p-orbitals. 

The sp 3 WFs of the N atom are not well suited as a 
starting point for applications of the Newns Anderson 
model, since they do not represent the energy levels of 
the free atom. This problem can be overcome by diago- 
nalizing the Hamiltonian matrix in the WF basis, within 
the subspace spanned by the four sp 3 orbitals. The result 
of the subspace diagonalization is a set of four atomic or- 
bitals consisting of one s-like and three p-like orbitals, 
each centered at the N atom. Two of the p-orbitals lie in 
the surface plane (the xy-plane) and are directed along 
the arrows shown in Fig. 1131 while the third is oriented 
along the surface normal (the z-axis). We shall refer to 
the p-orbitals as p x ,Py and p z , respectively. The ener- 
gies corresponding to the atomic orbitals are (in eV): 
s s = -14.8, e z = -2.4, e x = -3.7, and e y = -4.2. We 
notice, that the energy of the p x and p y orbitals differ 
even though the symmetry of the system suggests that 
they should be equal. The reason for this is that the 
WFs break the four-fold rotation symmetry of the sys- 
tem, i.e. the subspace spanned by the four sp 3 WFs 
is not invariant under the same symmetry transforma- 
tions as the Hamiltonian. This is not surprising, since 
the WFs are constructed solely from a criterion of max- 
imal localization and no attempts are made to conserve 
symmetries. On the other hand we have found that by 
increasing the parameter Eq above the value Eq = Ep 
used in the present example, the symmetry between p x 
and p y can be restored. The price one has to pay is that 
the copper s-like WFs become less localized due to the 
further constrains on the localization space implied by 
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the larger value of Eq. From the Hamiltonian in the WF 
basis we can also obtain the coupling, V, between each of 
the atomic nitrogen orbitals and its corresponding group 
orbital. These are quite similar and vary from 3.1 eV to 
3.8 eV. 

In Fig.^Jwe show the calculated PDOS for each of the 
three nitrogen p-orbitals (upper panel). Although the on- 
site energies of the p x and p y orbitals differ (as discussed 
above), their PDOS are rather similar. The PDOS of 
the corresponding group orbitals have been calculated 
with all coupling matrix elements to the N orbitals set 
to zero, i.e. the adsorbate states have effectively been 
decoupled from the surface. The result is shown in the 
lower panel of Fig. For all three orbitals, the on-site 
energies lie within the band. Due to the strong coupling, 
bonding and anti-bonding resonances are formed at the 
band edges around — 7 eV and eV as can be seen in 
the upper panel of the figure. This is the limit of strong 
chemisorptioniSl Since the four orbitals span all states 
with significant weight on the N atom, this representation 
provides a full representation of the nitrogen bonding. 



VI. CONCLUSIONS 

We have presented a practical method for construct- 
ing partly occupied WFs for a wide range of systems. 
The method employs a bonding-antibonding closing pro- 
cedure to filter out a set of unoccupied states, called the 
extra degrees of freedom, which serve to improve the lo- 
calization of the WFs. The determination of the extra de- 
grees of freedom is based on a minimization of the spread 
of the resulting WFs. We derived expressions for the gra- 
dients of the spread functional and showed how these can 
be combined with a Lagrange multiplier scheme to min- 
imize the spread functional. 

The generality of the scheme was demonstrated by ap- 
plying the method to three different systems. As an ex- 
ample of an isolated system, we considered a Sis cluster, 
and showed how different sets of WFs could be obtained 
by varying the number of extra degrees of freedom. A 
similar analysis was performed for a copper crystal, where 
we found results very similar to those of Souza et ali&. 
Finally we studied in detail the WFs of a Cu(100) surface 
with a nitrogen coverage of 0.5. In many cases we were 
able to obtain a special set of WFs with a particularly 
high degree of symmetry and localization, by maximizing 
the average spread of the WFs. Moreover, the condition 
of maximal average localization was shown to coincide 
with a complete matching of bonding and antibonding 
states. 



VII. ACKNOWLEDGMENTS 

We acknowledge support from the Danish Center for 
Scientific Computing through Grant No. HDW- 1101-05. 



APPENDIX A: SPREAD FUNCTIONAL FOR 
VANDERBILT ULTRASOFT 
PSEUDO-POTENTIALS 

For Vanderbilt ultra-soft pseudo-potentials 24 the op- 
timal smoothness of the pseudo-wavefunctions is ob- 
tained by relaxing the norm-conserving constrains for the 
pseudo-wavefunctions. This results in a generalized or- 
thonormality relation^ 

^ i \S\^ j ) = S ij . (Al) 

The Hermitian operator S is given by 

S = l + Y.Y.1™\P I n ){Pm I l (A2) 

/ nm 

where the index / denotes the atoms in the system, and 
q nm is given by 

q nm = J drQ£ m (r). (A3) 

The functions {/%} and {Q^ m } are all localized func- 
tions centered at atom /. The functions {Qn m } describe 
the augmentation charge not contained in the smooth 
pseudo-wavefunctions, and they must therefore be in- 
cluded in the calculation of the spread of the wavefunc- 
tions. 



1. Large supercells 

In the case of large supercells, using the T-point ap- 
proximation, Bernasconi and Madden^S derived the fol- 
lowing expression for the contribution to Za^ from the 
augmentation charges Q^ m (r): 

= £M/&K/£hM / dre- iG ^Q mn (r) (A4) 

2. Periodic systems 

For the periodic case, using a uniform k-point grid, we 
write Z™ s as 

Zl s = ^Z^ kk '. (A5) 

kk' 

Here again we use the exact correspondence between the 
Bloch states ip n k and the T-point eigenstates of the re- 
peated cell. In the repeated cell we use the notation, 

= W>ik|/#*> = (^k|/^ t= V k - Rt (A6) 

and 

QlL = Y.Q™{G)e- lG <*-^ (A7) 

G 
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Rt is here a real space translation vector, given in terms 
of the basis a, iiai + £ 2 a 2 + i3 a 3> t = (ti,*2>*3)- We 
will use h In = /i / ' t=0 ' n in what follows. Inserting hf^ 
and Qnm f rom Eqs. (|A6|) and (|A7|) . together with the 
Bloch states, ipiic, into the T-point expression for zi" s ' ^ 
in Eq. (|A4|) . we find 



kk ' = E E »SW / dre^-Q^M 

+■ r « 



t I.nm 



(A8) 



The sum over t is for t, = 0,..,Ni — 1, see Eq. (|2*Tll . 
Inserting the left hand side of Eqs. 1A7|1 and (|A6(I . and 

rearranging, we find 



7 (us,0),kk' 



E ^ikhj™! E e~ i ^ k_k '^' Rt E Qnm(G')e~ iG ' Rt / dre 
J,nm t g 



i(G-G a )-r 



(A9) 



Finally, we arrive at our expression for Z^'°''' kk 



7 («s,0),kk' 

J a,ij 



(A10) 



which is non-zero only when k and k' fulfills the condition the T-point formalism, Eq. I|A4(1 . as a special case, 
in Eq. I|29|) . Again we see that this expression contains 



1 G. H. Wannier Phys. Rev. 52, 191 (1937). 

2 N. W. Ashcroft and N. D. Mermin Solid State Physics 
(Saunders, New York, 1976). 

3 G. D. Mahan, Many-Particle Physics (Plenum Press, New 
York, 1990). 

4 S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999); C. K. Sky- 
laris et al., Phys. Rev. B 66, 035119 (2002); D. R. Bowler, 
T. Miyazaki and M. J. Gillan, J. Phys. Cond. Matt. 14, 
2781 (2002); J. L. Fattebert and J. Bernholc, Phys. Rev. 
B 62, 1713 (2000); J. Kim, F. Mauri and G. Galli, Phys. 
Rev. B 52, 1640 (1995). 

5 R. Resta Rev. Mod. Phys. 66, 899 (1994). 

6 D. Vanderbilt and R. D. King-Smith Phys. Rev. B 47, 1651 
(1993). 

7 A. Calzolari, N. Marzari, I. Souza, M. B. Nardelli Phys. 
Rev. B 69, 035108 (2004). 

8 K. S. Thygesen and K. W. Jacobsen, Chem. Phys. (in pro- 
duction), arXiv:cond- mat/0501238 

9 S. F. Boys, Rev. Mod. Phys. 32, 296 (1960). 

10 J. M. Foster and S. F. Boys, Rev. Mod. Phys. 32, 300 
(1960). 

11 C. Edmiston and K. Ruedenberg, Rev. Mod. Phys. 35, 457 
(1963). 

12 J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989) 

13 G. Berghold, C. J. Mundy, A. H. Romero, J. Hutter and 
M. Parrinello, Phys. Rev. B 61, 10040 (2000). 

14 P. L. Silvestrelli, N. Marzari, D. Vanderbilt and M. Par- 



rinello, Solid State Commun. 107, 7 (1998). 

N. Marzari and D. Vanderbilt Phys. Rev. B 56, 12847 

(1997) . 

K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. 

Rev. Lett. 94, 026405 (2005). 

P. L. Silvestrelli Phys. Rev. B 59, 9703 (1999). 

M. Iannuzzi and M. Parrinello, Phys. Rev. B 66, 155209 

(2002). 

I. Souza, N. Marzari and D. Vanderbilt, Phys. Rev. B 65, 
035109 (2001). 

L. Bernasconi and P.A. Madden, Jour, of Molecular Struc- 
ture (Theochem) 544,49,(2001) 

R. Resta and S. Sorella, Phys. Rev. Lett. 82, 370 (1999). 
H. Jonsson, G. Mills and K. W. Jacobsen, in Classical 
and Quantum Dynamics in Condensed Phase Simulations, 
(World Scientific, Singapore, 1998) 

B. Hammer, L.B. Hansen, and J.K. N0rskov, Phys. Rev. 
B 59, 7413 (1999); S.R. Bahn and K.W. Jacobsen, Comp. 
Sci. Eng. 4, 56 (2002); The Dacapo code can be down- 
loaded at http://www.camp.dtu.dk/campos 
D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 
H. Ohnishi, Y. Kondo and K. Takayanagi, Nature 395, 780 

(1998) . 

K. Raghavachari and V. Logovinsky, Phys. Rev. Lett. 55, 
2853 (1985). 

D. M. Newns, Phys. Rev. 178, 1123 (1969). 



